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We investigate viscous and non-viscous flow in two-dimensional self-affine fracture joints through 
direct numerical simulations of the Navier-Stokes equations. As a novel hydrodynamic feature of 
this flow system, we find that the effective permeability at higher Reynolds number to cubic order, 
falls into two regimes as a function of the Hurst exponent h characterizing the fracture joints. For 
h > 1/2, we find a weak dependency whereas for h < 1/2, the dependency is strong. A similar 
behavior is found for the higher order coefficients. We also study the velocity fluctuations in space 
of a passive scalar. These are strongly correlated on smaller length scales, but decorrelates on larger 
scales. Moreover, the fluctuations on larger scale are insensitive to the value of the Reynolds number. 

PACS numbers: 47.53.+n, 47.56.+r, 47.60.-K, 83.50.-v 



The transport of oil in carbonate reservoirs is domi- 
nated by flow through internal fractures. As about 10- 
15% of the world's known reservoirs consist of carbonates 
and that these contain about 50% of the remaining oil, 
it is surprising that this problem has received so little 
attention as it has from not only the physics community, 
but also the engineering communities. It is not uniquely 
from a practical, economical point of view that this is 
a worthwhile problem to study. The progress made in 
fracture morphology over the last years, has uncovered a 
large number of basic questions to be answered, and the 
present problem is one of these. It has become increas- 
ingly clear that the morphology of brittle fractures fol- 
lows certain scaling laws first observed in the mid-eighties 
, demonstrating that such surfaces are self-affine. That 
is, by rescaling the in-plane length scale x by Xx requires 
the out-of-planc scale y to be rescaled by X h y for the 
fracture surface to remain statistically unchanged, with 
h being the roughness exponent. In the early nineties, 
it was suggested that h is universal and has a value of 
0.8 2}. There is today some evidence that other values 
of the Hurst exponent may occur depending on the ma- 
terial and how the fracture has been produced [sj, thus 
suggesting more than one universality class. 

For single-phase fluid flow in porous media and frac- 
tures, it is common to characterize the system in terms of 
Darcy's law 0,0, which assumes that a global index, the 
permeability k, relates the average fluid velocity V with 
the pressure drop AP across the system, V = —kAP/^,L. 
Here L is the length of the sample in the flow direction 
and fi is the viscosity of the fluid. However, in order to 



understand the interplay between porous structure and 
fluid flow, it is necessary to examine local aspects of the 
pore space morphology and relate them to the relevant 
mechanisms of momentum transfer (viscous and inertial 
forces). Flow in self-affine faults was to our knowledge 
first discussed by Roux et al. , where they considered a 
fracture fault consisting of two matching walls that have 
been moved along the fracture plane with respect to each 
other. If the in-plane movement is x, then - due to the 
self-affinity of the fault - the amplitude of the fault will 
be proportional to x h . It was claimed in that the 
permeability of a two-dimensional fracture fault would 
scale as x 2h , based on using the self-affinity to deter- 
mine the scaling properties of an effective channel width. 
However, numerical studies by Gutfraind et al. Q using 
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FIG. 1: (a) Contour plot of the local velocity magnitude along 
a typical realization of the self-affine rough channel (h — 0.8) 
subjected to low Reynolds conditions. Fluid is pushed from 
left to right. The colors ranging from blue to red correspond 
to low and high velocity magnitudes, respectively. The solid 
line in the center of the channel corresponds to the streamline 
that divides the flux into two equal parts. The zoom shown 
in (b) reveals the presence of high velocity spots in the flow 
field induced by the occurrence of high slopes in the channel 
geometry. 
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FIG. 2: Dependence of the hydraulic resistance G = 
~APw 2 /fiVL on the Reynolds number for different values 
of the roughness exponent h. In all cases, the plateau corre- 
sponding to Darcy's law (constant G) is followed by a nonlin- 
ear regime that reflects the effect of convection on the flow. 
The error bars are smaller than the symbols and the solid 
lines are the best fit to the data of Eq. @ ■ 



FIG. 3: (a) Variation of the parameter a as in Eq. 1^1 with 
the exponent h characterizing the roughness of the channel ge- 
ometry. The two dashed lines correspond to different regimes 
dictated by the nature of the correlations. The crossover at 
h ~ 0.5 delimits anti-correlated (h < 0.5) and correlated 
(h > 0.5) signals. The inset shows that the parameters (5 
and 7 of Eq. @ behave in a very similar way. 



lattice gas methods, demonstrated that the permeability 
should be controlled by the narrowest neck in the system. 
In a subsequent work ||, the permeability of self-affinc 
rough fractures with wide and narrow apertures has been 
investigated analytically and confirmed through numeri- 
cal simulations with the lattice Boltzmann method. The 
flow in self-affine fracture joints, where the two matching 
fractures are moved apart by a distance w but without 
any relative movement in fracture plane (x = 0), has 
been previously studied by Skjctnc et al. 

Strictly speaking, the concept of permeability as a 
global index for flow in fractures should be restricted to 
viscous flow (linear) conditions. More precisely, Darcy's 
law should only be applicable for flow at sufficiently low 
Reynolds number, defined here as Re = pVw/fi, where 
w is the fracture opening. It is well known, however, 
that the role of inertial forces (convection) to flow in dis- 
ordered media should be examined in the framework of 
the laminar flow regime, before assuming that fully de- 
veloped turbulence effects are already present [j, [jjj, UM ■ 
Much less effort has been dedicated to address this prob- 
lem in the specific context of flow through rough frac- 
tures. For instance, the computational fluid dynamics 
simulations presented in Ref. [9j indicate that vanishing 
weak and strong inertial flows can be quantitatively de- 
scribed in different ways, namely the Darcy, weak inertia 
and Forchheimer empirical equations, respectively. Here 
we investigate by direct simulation of the Navier-Stokcs 
equations the departure from Darcy's law in laminar flow 
through self-affine fracture joints. We confirm that the 
physical description underlying a classical cubic equation 
provides a legitimate correlation for the flow in the frac- 
ture over a wide range of Reynolds number conditions. 
We then demonstrate that it is also possible to charac- 
terize this transition from linear to nonlinear behavior 
in terms of the spatial correlations in the fluid velocity. 



This allows us to elucidate certain characteristics of the 
fluid flow phenomenon in fracture fault geometries that 
have not been studied before. 

The self-affine surfaces are generated here with a 
Fourier method 0| f° r which we specify the length N 
and width m, both in terms of the number of nodes, 
and the roughness exponent h. The surface is given 
by yi (measured in units of the lattice constant), and 
m = max?/; — miny^. Since the lattice constant is 8, 
the length of the system is L = N8 and its amplitude is 
given by a = mS. When L is to be changed while 8 is 
to be kept fixed, we scale N — > XN and m — > Xin. On 
the other hand, when 8 is to be changed, while L is to 
be kept fixed, we let N — > XN and keep m fixed. The 
fracture opening w is kept fixed. 

The mathematical description for the fluid mechan- 
ics in the interstitial porous space is based on the as- 
sumptions that we have a continuum, Newtonian and in- 
compressible fluid flowing under steady state conditions. 
Thus, the Navier-Stokes and continuity equations reduce 
to 

/)u-Vu=-Vp + /iV 2 u, (1) 



V • u = , (2) 

where u and p are the local velocity and pressure fields, 
respectively, and p is the density of the fluid. In our sim- 
ulations, we consider non-slip boundary conditions along 
the entire solid-fluid interface. In addition, the changes 
in velocity rates are assumed to be zero at the exit x = L 
(gradientless boundary conditions), whereas a uniform 
velocity profile, u x (0,y) = V and u y (0,y) = 0, is im- 
posed at the inlet of the channel. 

The numerical solution of Eqs. Q and J5J) for the ve- 
locity and pressure fields in the pore space is obtained 
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through discretization by means of the control volume 
finite-difference technique |l3j | . Considering the complex 
geometries involved, we build an unstructured mesh of 
triangular grid elements based on a Delaunay network. 
For the largest fracture investigated, namely L = 512, 
around 10 6 cells are necessary to generate satisfactory 
results when compared with numerical meshes of larger 
resolution. The convergence criteria used in the simula- 
tions are defined in terms of residuals, i.e., the degree up 
to which the conservation equations are satisfied through- 
out the flow field. In all our simulations, convergence is 
considered to be achieved only when each of the residuals 
falls below 1CT 6 . 

In Fig. ^ we show the contour plot of the velocity mag- 
nitude in a typical self-affine channel under viscous flow 
conditions, i.e., at very low Reynolds number (Re = 10). 
Clearly, the spots of high velocity correspond to those re- 
gions with high slopes in the channel due to their reduced 
effective areas for flow, namely the cross-sections orthog- 
onal to the walls. Also shown in Fig. ^ is the streamline 
that divides the flow exactly into two zones of equal flux. 
As discussed later, this line will be adopted here as a 
reference location to study the spatial correlations in the 
local fluid velocity for different Reynolds conditions. 

The approach we adopt here to macroscopically char- 
acterize the effect of convection on flow through the self- 
affine channel is to employ the cubic relation 

AP , ■yp 2 V 3 
— = anV + f3pV 2 + — , (3) 

where the coefficient a corresponds to the reciprocal of 
the permeability of the porous material and the two re- 
maining terms containing f3 and 7 can be interpreted, 
respectively, as second and third order corrections that 
should account for the contribution of inertial forces in 
the fluid flow. At sufficiently low Reynolds, Eq. re- 
duces to Darcy's law. Rewriting {3} in terms of the 
Reynolds number we obtain, 

G = aw 2 + PwRe + 1 Re 2 , (4) 

with G = —APw 2 /p J VL being a dimcnsionless measure 
of the hydraulic resistance. Fig.|5]shows the results of our 
flow simulations in terms of the variables G and Re for 
different values of the roughness exponent h. After com- 
puting and averaging G over a total of 10 realizations for 
each value of h and a wide range of Reynolds numbers, 
we fit the results with Eq. (0} to estimate the coefficients 
a, (3 and 7. In agreement with real flow experiments, 
we observe a transition from constant G (Darcy's law) 
to nonlinear flow behavior at a value of Re that depends 
significantly on the roughness. As shown in Fig. |3 the 
coefficient a generally decreases with the roughness expo- 
nent h. There is, however, a clear crossover at a value of 
h w 0.5 separating two distinct regimes that characterize 
the influence of the channel geometry on the parameter 
a. Interestingly, the value h = 0.5 corresponds exactly 
to the transition point between anti-correlated (h < 0.5) 
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FIG. 4: (a) Profile of the velocity magnitude along the stream- 
line located at the center of the rough channel (see Fig. for 
a flow field calculated at low Reynolds conditions (Re = 10). 
(b) The same as in (a), but for Re = 1500. 

and correlated (h > 0.5) self-affine interfaces 01- To the 
best of our knowledge, this is the first time that a clear 
connection can be drawn in order to relate the degree of 
correlation in the interface geometry and the single-phase 
flow behavior through self-affine fractures. The results 
shown in the inset of Fig. [3] reveal that the coefficients (3 
and 7 also decrease with h and display similar crossovers 
between correlated and anti-correlated geometries. 

Next we study the fluctuations in velocity of a mass- 
less particle (i.e., a passive tracer) released right at the 
center of the channel inlet. As mentioned earlier, this 
particle will then follow a trajectory that coincides with 
the streamline dividing the flow into two regions of equal 
flux (see Fig. QJ. In Fig. we show the variation of 
its normalized velocity magnitude u* = u/V, along the 
main flow direction x in a typical realization of the rough 
channel and for two different values of the Reynolds num- 
ber. At low Re (Fig. U^), the location and intensity of 
the velocity peaks essentially correspond to the spatial 
variation in amplitude of the local slopes along the frac- 
ture. At high Re (Fig. ^Jd), the situation becomes quite 
different. Due to inertia, the effect on the flow field of 
the local channel geometry reveals a persistent behav- 
ior in the local velocity fluctuations, when compared to 
the results obtained at low Re (Fig. 0Ji). More precisely, 
whenever a sudden increase in velocity is observed due 
to the presence of a narrow constriction in the channel, 
the signal tends to decay more slowly at higher Re condi- 
tions, before the particle experiences another substantial 
amplitude fluctuation. It is interesting to note that the 
same sequence of peaks and valleys presented in Fig. [3Ji 
can also be observed in Fig. \%}p, but with the difference 
that the background velocity is generally much higher at 
high Re values because of inertial effects. 

Due to the interplay between flow and the self-affine 
characteristic of the fracture interfaces utilized here [l4| , 
one should expect the velocity profiles shown in Fig. 0] to 
contain a certain degree of correlation. More precisely, 
long-range power-law correlations in velocity magnitude 
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FIG. 5: (a) Logarithmic plot of the detrended fluctuation 
function F(Ax) calculated from 10 realizations of velocity 
magnitude profiles (see Fig. |1J and four different values of 
the Reynolds number. The two straight lines show the best 
power law fits to the data in the scaling regions. The dif- 
ference in the slopes indicates the passage from a highly cor- 
related (£ ~ 1.7) to an uncorrelated ~ 0.5) series. The 
inset shows that the function ^(Aa;) is invariant with Re for 
Re < 100. 

are identified and quantified here b y m eans of the de- 
trended fluctuation analysis (DFA) [15|. According to 
this method, one can avoid the spurious observation of 
long-range correlations due to non-stationarities by inte- 
grating the time series and mapping it to a self-affinc 
stochastic process. The velocity profile u*(x) is inte- 
grated, after subtracting the overall velocity average, and 
divided into boxes of equal length Ax. The local linear 
trend is then calculated as the least-squares straight line 
that fits the data within each interval. The difference 
between the integrated time series and the local trend 
in each box gives the detrended time series, from which 



we can compute the root-mean-squarc fluctuation F(Ax) 
|l5| . Scaling is present if F(Ax) has a power-law depen- 
dence on Ax, F(Ax) ~ (Ax) 1 * , where the exponent £ 
characterizes the nature of the long-range correlations. 

The results in Fig. [3] show that, regardless of the 
Reynolds number, the function F(Ax) displays a highly 
correlated power-law regime (£ = 1.7 ± 0.03) at small 
length scales followed by a typically uncorrelated scaling 
(C = 0.5±0.01)at larger values of the Ax window. The 
only difference is that, above a sufficiently high value of 
Re (see the inset of Fig.|SJ), the crossover from correlated 
to uncorrelated behavior starts to increase with Re. This 
is compatible with our previous qualitative analysis based 
on the simple inspection of the profiles in Fig. 0] The 
first (highly correlated) regime is a direct consequence of 
the channel self-affinity inducing long-range correlations 
in tracer velocity, where a large exponent denotes the 
fractional-Brownian aspect of the signal. 

In summary, we have presented a numerical study of 
single-phase flow in self-affine fracture joints. We find 
that each of the coefficients a, (3 and 7, characterizing 
the cubic generalization of the Darcy equation J2J de- 
pends nearly linearly on the Hurst exponent h in two 
distinct regimes that meet at h « 0.5. We also follow 
a passive scalar particle released at the center of the 
channel, finding that its velocity fluctuations are char- 
acterized by strong spatial correlations on small scales, 
crossing over to uncorrelated random-walk-like behavior 
on larger scale. The large-scale behavior is insensitive to 
Re. 
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